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Abstract 

It  is  known  that  small  relative  perturbations  in  the  entries  of  a  bidiagonal  matrix 
only  cause  small  relative  perturbations  in  its  singular  values,  independent  of  the  values 
of  the  matrix  entries.  In  this  paper  we  show  that  a  matrix  has  this  property  if  and  only 
if  its  associated  bipartite  graph  is  acyclic.  We  also  show  how  to  compute  the  singular 
values  of  such  a  matrix  to  high  relative  accuracy.  The  same  algorithm  can  compute 
eigenvalues  of  symmetric  acyclic  matrices  with  tiny  componentwise  relative  backward 
error.  This  class  includes  tridiagonal  matrices,  arrow  matrices,  and  exponentially  many 
others. 

1     Introduction 

In  [9]  it  was  shown  that  small  relative  perturbations  in  the  entries  of  a  bidiagonal  matrix 
B  only  cause  small  relative  perturbations  in  its  singular  values.  This  is  true  independent 
of  the  values  of  the  nonzero  entries  of  B.  This  property  justifies  trying  to  compute  the 
singular  values  of  B  to  high  relative  accuracy,  and  is  essential  to  the  error  analyses  of  the 
corresponding  algorithms  [9]. 

Since  this  attractive  property  of  bidiagonal  matrices  is  independent  of  the  values  of  the 
nonzero  entries,  it  is  really  just  a  function  of  the  sparsity  pattern  of  bidiagonal  matrices. 
In  this  paper  we  completely  characterize  those  sparsity  patterns  with  the  property  that 
independent  of  the  values  of  the  nonzero  entries,  small  relative  perturbations  of  the  matrix 
entries  only  cause  small  relative  perturbations  of  the  singular  values.  The  characterization 
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is  simple:  a  sparsity  pattern  has  this  property  if  and  only  if  its  associated  bipartite  graph 
is  acyclic. 

We  define  this  graph  as  follows.  Let  5  be  a  sparsity  pattern  for  m  by  n  matrices;  in 
other  words,  S  is  a  list  of  the  entries  permitted  to  be  nonzero.  Let  G(S)  be  a  bipartite  graph 
with  one  group  of  nodes  {rl5...,rm}  representing  the  m  rows  and  one  group  {ci,...,cn} 
representing  the  n  columns.  There  is  an  edge  from  r,-  to  Cj  if  and  only  if  AtJ-  is  permitted 
to  be  nonzero.  (We  wiD  sometimes  write  G(A)  instead  of  G(S),  where  S  is  the  sparsity 
pattern  of  A.) 

We  also  present  another  perturbation  property  of  acyclic  matrices  which  is  quite  strong: 
multiplying  any  single  matrix  entry  by  any  factor  f3  ^  0  cannot  change  any  singular  value 
by  more  than  a  factor  of  j3  (either  up  or  down). 

Sparsity  patterns  with  this  property  have  at  most  n  +  m  —  1  nonzero  entries.  There 
art'  a  great  many  such  sparsity  patterns.  Let  us  consider  only  m  by  n  sparsity  patterns  S 
which  cannot  be  permuted  into  block  diagonal  form  (this  means  G(S)  is  connected).  Then 
the  number  of  different  such  sparsity  patterns  is  equal  to  the  number  of  spanning  trees  on 
connected  bipartite  graphs  with  m  +  n  vertices;  this  number  is  m^^n™-1  [5,  p.  38]  [3].  If 
we  only  wish  to  count  sparsity  patterns  which  cannot  be  made  identical  by  reordering  the 
rows  and  columns,  a  very  simple  lower  bound  on  the  number  of  such  equivalence  classes  is 
m"~17?m-1/(n!m!).  In  the  square  case  n  =  m,  Stirling's  formula  lets  us  approximate  this 
lower  bound  by  e2n /(2i:n3),  which  grows  quickly. 

Since  we  know  the  singular  values  of  these  acyclic  matrices  are  determined  to  high 
relative  accuracy  by  the  data,  it  makes  sense  to  try  to  compute  them  this  accurately. 
We  present  a  bisection  algorithm  which  does  this.  The  same  algorithm  can  compute  the 
eigenvalues  of  arbitrary  "symmetric  acyclic"  matrices  with  tiny  componentwise  relative 
accuracy.  We  define  symmetric  acyclicity  of  a  symmetric  matrix  as  follows.  Given  a  sparsity 
pattern  S  of  an  n  by  n  symmetric  matrix,  we  define  a  graph  G'(S)  by  taking  n  nodes,  and 
connecting  node  i  to  node  j  ^  i  if  and  only  if  the  (i,j)  entry  is  nonzero.  The  symmetric 
sparsity  pattern  S  is  called  "symmetric  acyclic"  if  the  graph  G'(S)  is  acyclic.  (We  will 
sometimes  write  G'(A)  instead  of  G'(S)  where  S  is  the  sparsity  pattern  of  A.)  The  algorithm 
evaluates  the  inertia  of  such  a  matrix  by  doing  symmetric  Gaussian  elimination,  with  the 
older  of  elimination  determined  by  a  postorder  traversal  of  G'(S). 

In  summary,  the  well-known  attractive  properties  of  bidiagonal  matrices  B  and  symmet- 
ric tridiagonal  matrices  T,  that  the  singular  values  of  B  can  be  computed  to  high  relative 
accuracy  and  the  eigenvalues  of  T  computed  with  tiny  componentwise  relative  backward  er- 
ror, have  been  extended  to  "acyclic"  matrices.  In  the  case  of  computing  singular  values,  we 
have  shown  that  this  extension  is  complete:  no  other  sparsity  patterns  have  this  property. 
We  strongly  suspect  that  the  set  of  symmetric  acyclic  matrices  is  also  the  complete  set  of 
symmetric  matrices  whose  eigenvalues  can  be  computed  with  tiny  componentwise  relative 
backward  error  independent  of  the  values  of  the  matrix  entries. 

Other  algorithms  for  the  special  case  of  "arrow"  matrices  are  discussed  in  [1,2,15,22]. 
This  work  generalizes  the  adaptations  of  bisection  to  arrow  matrices,  and  is  almost  certainly 
more  stable  than  the  QR  based  schemes. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  states  the  perturbation  theorem 
for  the  singular  values  of  acyclic  matrices,  and  section  3  proves  it.  Section  4  shows  how 
to  compute  eigenvalues  of  symmetric  acyclic  matrices  with  tiny  componentwise  relative 


backward  error,  and  applies  this  to  compute  the  singular  values  of  acyclic  matrices  to  high 
relative  accuracy.  Section  5  give  some  examples  of  matrices  with  acyclic  sparsity  patterns. 
Section  6  discusses  algorithms  and  open  problems. 

2      Statement  of  Perturbation  Theorem  for  Singular  Values 

In  this  section  we  define  three  properties  of  sparsity  patterns  of  matrices,  one  about  graph 
theory  and  two  about  perturbation  theory.  Our  main  result,  which  we  prove  in  the  next 
section,  is  that  these  properties  are  equivalent. 

Le1  A  be  an  m  by  n  matrix  with  a  fixed  sparsity  pattern  S. 

Property  1    .  G(S)  is  acyclic. 

Property  2  .  Given  sparsity  pattern  S ,  there  exist  ]x>sitive  constants  eo  and  (  with  the 
following  property.  Let  A  be  any  matrix  with  sparsity  pattern  S,  and  A{:  any  nonzero  entry. 
( 'hoosi  any  |e|  <  to,  and  let  A'  =  A  except  for  A\;  =  Aij(l  +  e).  Then  for  all  singular  values 
trk(A') 

(i  -  CklK(^)  <  °k(A')  <  (i  +  CklK(^) 

///  other  words  a  sufficiently  small  relative  perturbation  e  in  any  single  matrix  entry  cannot 
caust   a  relativt  pi  rturbation  greater  than  (t  in  any  singular  value. 

Up  entries  of  A  are  simultaneously  perturbed,  Property  2  can  be  applied  p  times  to  show 
no  singular  value  can  change  by  afactor  outside  the  interval  from  (1— Ckl)p  =  l—p(\€\  —  0(e2) 
to  ( J  +  (,'|c| )''  =  1  +p(|e|  -  0{e2).  Property  2  seems  rather  weak,  since  it  imposes  no  bounds 
on  £q  1io>  C-  Still,  since  eo  and  £  are  independent  of  the  matrix  entries,  it  is  actually 
demanding  quite  a  bit  of  S.  The  last  property  is  even  stronger: 

Property  3  .  Given  sparsity  pattern  S,  let  A  be  any  matrix  with  this  sparsity,  and  Aij  any 
nonzero  i  ntry.  Let  (3  be  any  nonzero  constant.  Let  A'  —  A  except  for  A'{j  =  fiAir  Then  for 
all  singular  values  o~k(A') 

mto(|0|,|j8-1|)<7*(ii)  <  ak{A')  <  max(|/?|,  \firl\)*k(A) 

Property  3  is  much  stronger  than  Property  2  because  it  imposes  no  limit  to  on  the  size 
of  the  relative  perturbation,  and  because  it  asserts  (  =  1,  i.e.  that  the  relative  change  in 
the  singular  values  cannot  exceed  the  relative  change  in  the  single  perturbed  matrix  entry. 
In  the  case  of  simultaneous  small  relative  perturbations  of  size  at  most  /3  =  1  +  e  in  p  entries 
of  .4.  Property  3  implies  that  no  singular  value  can  change  by  a  factor  outside  the  interval 
from  ( 1  -  |e|  )p  =  1  -  p\e\  +  (e2)  to  (1  -  |£|)~p  =  1  +  p\e\  +  (e2).  Since  the  maximum  number 
of  nonzeros  is  m  +  n  —  1,  this  relative  perturbation  is  bounded  by  (m  +  n  —  l)|e|  +  0(e2). 

Our  main  result  is 

Theorem  1   Properties  1,  2  and  3  of  a  sparsity  pattern  S  are  equivalent. 


Figure  1:  Computing  DT  and  Dc 


if  q  is  a  row  node  then 
suppose  q  =  r, 
if  i\  is  the  root  then 

Dr,u  =  1 
else 

suppose  Cj  is  the  parent  of  q 

L>r,u    =    f/(  AijlJcJj  ) 

end 
else  [q  must  be  a  column  node)  then 
suppose  q  =  c} 
if  c t  is  the  root  then 

'  DeJj  =  1 
else 

suppose  r,  is  the  parent  of  q 
DcJj=  l/(AijDriii) 
end 
end  il 


3      Proof  of  Perturbation  Theorem  for  Singular  Values 

The  proof  of  equivalence  will  consist  of  the  following  steps.  We  already  know  that  Property  3 
implies  Property  2,  so  it  will  suffice  to  prove  Property  1  implies  Property  3,  and  Property  2 
implies  Property  1. 

Lemma  1  Let  A  have  sparsity  pattern  S,  and  suppose  G(S)  is  acyclic.  Then  there  are  di- 
agonal matrices  Dr  and  Dc  such  that  each  entry  of  DrADc  is  either  0  or  1.  Each  diagonal 
<  ntry  of  D,  or  Dc  is  a  quotient  of  monomials  in  the  entries  of  A.  In  each  monomial  each 
distinct  factor  AX}  which  appears  has  unit  exponent.  Each  A{j  can  appear  only  in  in  numer- 
ators of  entries  of  DT  and  denominators  of  entries  of  Dc,  or  vice  versa,  in  denominators 
of  entries  of  Dr  and  numerators  of  entries  of  Dc. 

Proof  Since  G(S)  is  acyclic,  it  is  a  forest  of  trees.  We  may  consider  each  tree  indepen- 
dently. We  traverse  each  tree  via  depth  first  search,  and  execute  the  program  in  Figure  1 
when  first  visiting  node  q. 

The  depth  first  search  visits  each  node  once.  Since  the  graph  is  bipartite,  row  nodes  and 
column  nodes  alternate,  so  the  parent  of  a  row  node  is  a  column  node  and  vice  versa.  Since 
each  node  is  visited  once,  the  above  program  is  executed  once  for  each  edge  in  the  tree,  i.e. 
once  for  each  nonzero  entry  A{j,  corresponding  to  the  edge  connecting  nodes  r\  and  Cj.  Thus 
each  Dr,u  and  Dc,jj  is  set  exactly  once.  Since  the  i,j  entry  of  DTADC  is  Dr,iiAijDc,jj,  we 
see  immediately  from  the  way  Dr<ll  and  Dcjj  are  defined  that  this  quantity  is  1  if  A{j  ^  0 
land  0  otherwise).   Since  each  Aij  is  used  once  during  the  graph  traversal,  each  DT>a  and 


I)  ti  iinisi  be  be  a  quotient  of  monomials.  If  AtJ  is  first  used  in  DT<a,  then  the  formulas  in 
i  he  above  program  and  the  fact  the  row  and  column  nodes  alternate  mean  that  A\3  will  only 
appear  in  denominators  of  entries  of  Dr  and  numerators  of  entries  of  Dc.  Alternatively,  if 
A,j  is  first  used  in  Dcjj,  then  A{j  will  only  appear  in  denominators  of  entries  of  Dc  and 
numerators  of  entries  of  Dr.  □ 

The  rest  of  the  proof  mimics  that  of  [4,  Thm.  1].  Let  E  be  the  matrix  of  ones  and  zeros 
with  sparsity  S,  so  that  DrADc  =  E.  Write  DT  =  Sr\Dr\  where  \Dr\  is  the  matrix  of  absolute 
values  of  D, .  and  ST  is  a  diagonal  matrix  with  \Sr\  =  I.  Similarly  write  Dc  =  SC\DC\.  Then 


A  =  D 


-i 


ed;1  =  sr 


-1 


\Dr\-lE\Dc\-lS;1  =  Sr 


-l 


\A\S, 


-l 

c 


so  that  .4  is  related  to  \A\  be  pre-  and  postmultiplication  by  diagonal  orthogonal  matrices. 
In  particular,  A  and  \A\  have  the  same  singular  values.  We  will  henceforth  assume  without 
loss  of  generality  that  A  is  nonnegative  and  so  DT  and  Dc  are  also  nonnegative. 

It  is  known  that  the  singular  values  of  A  are  the  same  as  the  positive  eigenvalues  of  the 
pencil 

0      A 
AT     0 


XI 


which  are  in  turn  the  same  as  the  positive  eigenvalues  of  the  equivalent  symmetric  definite 


pencil 

0 


0 


0      A 

AT     0 


A-/ 


Dr      0 
0      Dc 


0      E 
ET     0 


-A 


D2      0 
0      D2 


=  F-XD2 


Now  suppose  we  perturb  A  by  changing  nonzero  entry  AX]  to  0A{j,  resulting  in  the 
perturbed  matrix  .4'.  Apply  the  algorithm  in  Lemma  1  to  compute  a  new  D'r  and  D'c.  Since 
by  Lemma  1  the  entries  of  D'r  and  D'c  are  quotients  of  monomials  where  each  independent 
factor  appears  at  most  once,  each  entry  D'Tkk  must  equal  either  Dr,kk,  PDr<kk  or  fi~lDr,kk- 
An  analogous  statement  about  D'c  kk  and  Dc^  is  true.  Since  a  factor  Atj  must  appear  either 
in  numerators  of  Dr  and  denominators  of  Dc,  or  in  denominators  of  Dc  and  numerators  of 
D, .  we  have  two  cases: 

I.  Either  D'rkk  =  DrM  or  D'tM  =  $DrM,  and  either  D'ckk  =  DcM  or  D'ckk  =  f3~lDcM. 

■2.  Either  D'rkk  =  Dr>kk  or  D'T<kk  =  (3~1DtM,  and  either  D'ckk  =  DcM  or  D'ckk  =  0DCtkk. 

Note  we  may  multiply  Dr  by  any  nonzero  7  and  divide  Dc  by  7  without  changing  the 
fact  that  D,ADC  =  E.  Corresponding  to  the  above  two  cases,  we 

1.  divide  DT  by  \(3\ll2  and  multiply  Dc  by  I/?!1/2,  or 

2.  divide  Dc  by  \l3\^2  and  multiply  Dr  by  \p\ll2. 

The  end  result  will  be  D'r  and  D'c  matrices  each  of  whose  entries  differs  from  the  corre- 
sponding entry  of  Dr  and  Dc  by  factors  of  l/?^1/2.  In  particular,  this  implies 

I/3I"1  <  ^T^-,  <  \fi\    and     \f3\~'  <  4^  <  \p\ 


x1  D'ix 


xTD'2cx 


for  any  nonzero  vector  x.  Let  D'  =  diag  (D'lyD2)  as  we  above  defined  D  —  diag  (Di,D2). 
Then 

1-1 


yTD2y 


and 


for  any  nonzero  vector  y.  We  may  now  apply  [4,  Lemma  2]  to  conclude  that 

ak(A)  =  rain       m^     ^^~ 

X     J^X 

a*(A')  =  min       max  2      , 

S*      x  e  S*      x   D    x 

INl2  =  l 

where  the  minima  are  over  all  k  +  max(n,m)  dimensional  subspaces  Sk,  can  differ  by  no 
more  than  a  factor  of  f3.  This  proves  that  Property  1  implies  Property  3. 

Lemma  2  Let  A  have  sparsity  pattern  S,  and  let  all  its  nonzero  entries  be  independent 
indeterminates.  Then  G(S)  is  acyclic  if  and  only  if  all  minors  of  A  are  either  0  or  mono- 
mials. 

Proof  We  begin  by  noting  that  to  each  term  in  the  determinant  of  an  5  by  5  square 
matrix  M  corresponds  a  unique  perfect  matching  in  graph  G(M).  This  is  because  each 
term  in  t lie  determinant  corresponds  to  a  choice  of  s  entries  of  M  located  in  disjoint  rows 
and  columns,  and  each  such  choice  of  s  entries  selects  a  perfect  match  in  G(M). 

Now  suppose  a  square  submatrix  M  of  A  has  at  least  two  terms  in  its  determinant. 
These  correspond  to  two  different  perfect  matchings.  Take  the  symmetric  difference  of  the 
edges  in  these  matchings.  This  symmetric  difference  forms  a  cycle,  which  we  get  by  following 
edges  of  the  two  matchings  in  alternation.  Thus  G(M)  contains  a  cycle,  and  so  must  G{A) 
since  it  includes  G(M). 

Now  suppose  G(A)  contains  a  cycle.  Assume  without  loss  of  generality  that  it  is  a  simple 
cycle,  i.e.  it  is  connected  and  visits  each  node  once.  Let  M  by  the  corresponding  square 
submatrix.  This  cycle  determines  two  perfect  matchings  in  G(M),  consisting  of  alternate 
edges  of  the  cycle.  This  means  det(A/)  has  at  least  two  terms.  D 

To  prove  that  Property  2  implies  Property  1,  we  will  show  the  contrapositive.  So  assume 
G(A)  contains  a  cycle,  and  let  M  be  an  s  by  s  submatrix  whose  determinant  has  at  least 
2  terms.  This  means  we  may  choose  all  the  entries  of  M  to  be  nonzero  but  such  that 
M  is  exactly  singular.  Thus  its  singular  values  include  at  least  one  which  is  exactly  zero. 
Scale  M  so  that  its  entry  of  smallest  absolute  value  is  1,  and  let  a  =  \\M\\2  >  1.  Now  let 
A(AJ . )])  denote  the  matrix  with  sparsity  S,  submatrix  M,  and  other  nonzero  entries  equal 
to  //.  Then  A(M,0)  will  have  at  least  min(m,  n)  —  s  +  1  zero  singular  values,  min(m,n)  -  s 
from  the  zero  rows  and  columns  outside  M,  and  1  from  the  singularity  of  M.  By  standard 
perturbation  theory  A(M, n)  will  have  at  least  min(m,n)  -5+1  singular  values  no  larger 
t  ban  mni].  Now  change  a  smallest  entry  of  M  from  1  to  1  +  x  to  get  Mx;  thus  x  is  also  the 


relative  change  in  this  entry.  Then  |det(Mr)|  >  x,  and  so  <7mtn(MT)  >  \x\/(o  +  x)s  1.  This 
means  <ra(A(Mx,  T)))  >  \x\/(o  +  x)s+1  -  mnr],  whereas  as(A(M,  n))  <  mnr).  Thus 

as(A(Mx,r}))  >   {a4yn  ~  mnrj  x 

<rs(/l(.M,7?))    ~~  mnT/  mnn(o  +  x)s+1 

11  Property   1   held,  then  we  would  be  able  to  find  en  >  0  and  (  >  0  such  that  for  all 
0  <  .r  <  (a  and  7/  >  0  the  following  inequality  would  hold: 

i<C  • 


mnr)(a  +  x)s+l 

Since  we  can  make  77  as  small  as  we  like,  this  inequality  cannot  hold  for  any  finite  £•  Thus 
Property  2  cannot  hold.  This  completes  the  proof  that  Property  2  implies  Property  1,  and 
so  also  completes  the  proof  of  Theorem  1. 

4      A   bisection   algorithm   for   computing   eigenvalues   with 
tiny  backward  error 

Let  £M  denote  the  machine  precision.  We  will  assume  the  usual  model  of  floating  point 
error,  //(</  b)  =  (a  Q  6)(1  +  6)  with  |<5|  <  eM,  and  assume  neither  underflow  nor  overflow 
occur.  (Of  course,  a  practical  algorithm  would  need  to  account  for  overflow.  This  can 
he  done  analogously  to  the  way  overflow  is  accounted  for  in  standard  tridiagonal  bisection 
[13].) 

In  this  section  we  will  show  how  to  compute  the  eigenvalues  of  a  symmetric  acyclic 
matrix  T  with  tiny  componentwise  relative  backward  error.  Our  main  result  is 

Theorem  2  The  algorithm  in  Figure  2  computes  count(T,  x),  the  number  0}  eigenvalues  of 
l'  less  than  x,  with  a  backward  error  6T  with  the  following  properties: 

\6Tij\  <  (1.5w  +  2.5)£M|TtJ|  when  i  f  j. 

|6r,-,-|<(2t;  +  2)£M|x|.  I 

Hen  v  <  n  —  1  is  the  maximum  degree  of  any  node  in  the  graph  ofT.  In  other  words,  the 
computed  count(T.  x)  is  the  exact  value  o/count(T  +  6T,x)  where  ST  is  bounded  as  above. 

This  is  essentially  identical  to  the  standard  error  analysis  of  Sturm  sequence  evaluation 
for  symmetric  tridiagonal  matrices  [9,  Sec.  6]  [13]  (this  is  stronger  than  the  result  in  [20,  p. 
303]). 

Our  algorithm  simply  performs  symmetric  Gaussian  elimination  on  T  —  xl:  P(T  — 
xI)P  -  LDLT  where  P  is  a  permutation  matrix,  L  is  unit  lower  triangular  and  D  is 
diagonal.  Then  count(T, x)  is  simply  the  number  of  negative  diagonal  entries  of  D,  by 
Sylvester's  Inertia  Theorem  [16].  The  order  of  elimination  is  the  same  as  a  postorder 
traversal  of  the  nodes  of  the  acyclic  graph.  Since  leaves,  which  have  degree  1,  are  eliminated 
first,  there  is  no  fill-in  during  the  elimination,  and  all  off-diagonal  entries  L{j  of  L  can  be 
computed  by  simply  dividing  LtJ  =  T{j/Djj. 


Figure  2:  Computing  counter,  2) 

call  ci\t(i,d,s,x)  where  i  is  any  node  1  <  i  <  n 
return  count(T,  x)  =  s 

procedure  cnt(?,rf, s,x) 

/*  i  and  x  are  input  parameters,  d  and  5  are  output  parameters  */ 

d  =  Tu  -  x 

s  =  0 

for  all  children  j  of  i  do 

call  cnt(j,  d',s',x) 

d=d-  T?/d' 

s  =  s  +  s' 
end  for 

if  d  <  0,  then  3  =  6  +  1 
return  d  and  5 
end  procedure 


We  assume  the  graph  G\S)  is  connected,  since  otherwise  the  matrix  can  be  reordered  to 
be  block  diagonal  (one  diagonal  block  per  connected  component  of  G'(S)),  and  the  inertia 
of  each  diagonal  block  can  be  computed  separately.  The  algorithm  cnt(i,  d,  s,x)  in  Figure  2 
assumes  the  matrix  is  stored  in  graph  form.  Subroutine  cnt(i,d,s,x)  does  a  postorder 
traversal  of  the  acyclic  graph  G'(S),  and  may  be  called  starting  at  any  node  1  <  i  <  n.  In 
addition  tu  i.  x  is  an  input  parameter.  The  variables  d  and  s  are  output  parameters;  on 
return  .^  is  the  desired  value  of  count (T, a;). 

To  prove  Theorem  2,  we  will  exploit  the  acyclicity  of  T  to  show  that  each  computed 
quantity  and  original  entry  of  T  is  used  (directly)  just  once  during  the  entire  computation, 
and  then  use  this  to  "push"  the  rounding  error  back  to  the  original  data. 

We  see  that  each  entry  of  T  is  used  just  once  as  follows.  Tu  is  only  used  when  visiting 
node  /',  and  T%3  is  used  only  once,  when  visiting  i  if  j  is  a  child  of  i  or  when  visiting  j  if  i  is 
a  child  of  j  in  the  postorder  traversal  tree. 

Now  denote  the  d  computed  when  visiting  node  i  by  d±.  The  floating  point  operations 
performed  while  visiting  node  i  are  then 


/ 


\ 


T\ 


(4.1) 


di  =  fl    Tu  -  x  -  Y^ 

all  children 
V  j  of  i  j 

To  analyze  this  formula,  we  will  let  subscripted  e's  denote  independent  quantities 
bounded  in  absolute  value  by  eM.  We  will  also  make  standard  approximations  like 
(I  +£i)±1(l  +  £2)±1  =  l  +  2£3. 


Since  we  do  not  know  the  number  of  terms  or  the  order  of  the  sum  in  equation  (4.1),  we 
will  make  the  worst  case  assumption  that  there  are  v  <  n—  1  terms  where  v  is  the  maximum 
degree  ol  any  node  in  the  graph  G'(S).  This  leads  to 

T2 
di  =  (1  +  (v  +  l)£ia)Tu  -  (1  +  (v  +  l)eib)x  -  J2         (1  +  (v  +  3)e<i)-^-        (4.2) 

all  children 
j  of  i 


or 


all  children 
j  of  i 

Let  ;,„  be  the  roundoff  error  corresponding  to  £,a  committed  when  computing  dy  Then 
i =  r,,-s  +  (2„  +  2)£,s-  £  '"^"'t"^^''    (4.4) 

l  +  (»+l)£i„  ^-  rf,/(l+(l.+  l)£ja) 

all  children 


j  of  i 


or.  finally 


¥      t  ^o    x^  Y-  (d  +  (1.5t;  +  2.5)gl-J-»)ri-j)2 

all  children 
j  of  i 

where  d\  =  </,-/(  1  +  fJQ).  Equation  (4.5)  tells  us  that  the  d[  are  the  exact  diagonal  entries  of 
D  in  P[T  +  ST  -  xI)PT  =  LDLT .  Since  they  obviously  have  the  same  signs  as  the  cf,-,  this 
proves  Theorem  2. 

The  proof  depends  strongly  on  there  not  being  any  fill-in  and  on  each  off  diagonal  entry 
being  computable  by  a  single  division.  Since  these  properties  hold  if  and  only  if  the  graph 
Ci'(T)  is  symmetric  acyclic,  we  strongly  suspect  that  this  is  the  only  class  of  matrices  whose 
eigenvalues  can  always  be  computed  with  tiny  componentwise  relative  backward  error. 

We  now  apply  Theorem  2  to  compute  singular  values  of  acyclic  matrices  to  high  relative 
accuracy.  So  suppose  £  is  a  matrix  whose  graph  G(B)  is  acyclic.  Consider  the  symmetric 
matrix 

0      B 
BT     0 


A  = 


It  is  well  known  that  the  positive  eigenvalues  of  A  are  the  singular  values  of  B.  It  is 
also  immediate  that  the  graph  G'(A)  =  G{B).  Therefore  B  is  acyclic  if  and  only  if  A  is 
symmetric  acyclic,  so  we  can  apply  the  above  algorithm  to  compute  all  f?'s  singular  values 
to  high  relative  accuracy. 

One  other  algorithm  is  worth  mentioning.  If  A  is  symmetric  positive  definite  and  sym- 
metric acyclic,  then  its  Cholesky  factor  L  is  acyclic,  has  the  "lower  half"  of  the  sparsity 


pattern  of  A.  and  may  be  computed  by  using  algorithm  cnt.  It  may  occasionally  be  more 
accurate  to  compute  A's  eigenvalues  by  first  computing  L,  computing  its  singular  values  by 
bisection,  and  then  squaring  the  singular  values  to  get  A's  eigenvalues  [4].  This  is  the  case, 
for  example,  for  the  tridiagonal  matrix  with  2's  on  the  diagonal  and  l's  on  the  off-diagonal. 


5      Examples 

We  give  various  examples  of  acyclic  sparsity  patterns,  beginning  with  acyclic  G(S).  Given 
any  acyclic  sparsity  pattern,  others  can  be  generated  either  by  permuting  rows  and/or 
columns,  or  by  adding  more  zeros.  Since  all  square  acyclic  matrices  have  monomial  (or 
zero)  determinants,  this  means  we  can  permute  them  to  be  upper  triangular.  In  addition 
to  bi  diagonal  matrices,  some  other  examples  are 


x 

x 

x  x 

x     x 

x 


and 


x 

x 

x 

x  x 

x    x 

x 


1  = 


lb  get    symmetric  acyclic   matrices  A,   one  can  always  take  an  acyclic  B  and  set 
i  -  XI .  Some  other  examples  are 


x 

X 

X 

X 

X 

X 

X 

and 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

6     Algorithms  and  Open  Problems 

In  [8]  a  perturbation  theorem  for  singular  vectors  of  bidiagonal  matrices  is  proven,  which 
shows  that  the  appropriate  condition  number  for  the  z'-th  singular  vector  is  the  reciprocal 
of  the  relative  difference  between  the  i-th  singular  value  and  next  closest  one.  It  would  be 
interesting  to  extend  this  to  the  acyclic  case. 

Given  the  perturbation  theory,  it  would  be  nice  to  compute  the  singular  vectors  as 
accurately  as  they  deserve.  A  natural  candidate  is  inverse  iteration,  but  even  in  the  simple 
case  of  symmetric  tridiagonal  matrices,  open  problems  remain.  In  particular  there  is  no 
absolute  guarantee  that  the  computed  eigenvectors  are  orthogonal,  although  in  practice  the 
algorithm  can  be  made  quite  robust  [11]. 

In  the  "extreme"  cases  of  tridiagonal  and  arrow  matrices,  we  know  how  to  compute  the 
inertia  in  O(logn)  time,  using  the  so-called  parallel-prefix  algorithm  in  the  tridiagonal  case 
[17.19]  .  and  more  simply  in  the  arrow  case.  The  stability  in  the  tridiagonal  case  is  unknown, 
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bill  in  practice  it  appears  to  be  stable.  We  can  extend  this  to  the  general  symmetric  acylic 
case  in  two  ways.  First,  the  tree  describing  the  expression  whose  final  value  is  d,  has  at  most 
//  leaves.  From  [6]  we  know  any  such  expression  tree  can  be  evaluated  in  at  most  4log2  n 
parallel  steps,  although  stability  may  be  lost.  Another  approach,  which  includes  parallel 
prefix  and  the  algorithm  in  [15]  as  special  cases,  is  based  on  [14].  The  idea  is  to  simply 
evaluate  the  tree  greedily,  summing  k  leaves  of  a  single  node  in  0(\og2k)  steps  whenever 
possible,  and  collapsing  a  chain  of  k  nodes  into  a  single  node  via  parallel  prefix  in  0(log2  k) 
steps  whenever  possible.  If  we  could  understand  the  numerical  stability  of  parallel  prefix, 
we  could  probably  analyze  this  more  general  scheme  as  well. 

Divide  and  conquer  [7,10,18,12]  has  been  widely  used  for  the  tridiagonal  eigenproblem 
and  bidiagonal  singular  value  decomposition.  This  can  be  straightforwardly  extended  to 
the  acyclic  case.  In  terms  of  the  tree,  just  remove  the  root  by  a  "rank  one  tearing",  solve 
the  independent  child  subtrees  recursively  and  in  parallel,  and  merge  the  results  by  solving 
the  secular  equation  [21].  Any  node  can  be  the  root,  and  to  be  efficient  it  is  important  that 
no  subtree  be  large.  In  the  tridiagonal  case,  there  are  always  two  subtrees  of  nearly  equal 
size.  In  a  general  tree  one  can  only  make  sure  that  no  subtree  has  more  than  half  the  nodes 
of  the  original  tree  (this  is  easily  done  in  O(n)  time  via  depth  first  search). 

QR  does  not  appear  to  extend  beyond  the  tridiagonal  case.  The  case  of  arrow  matrices 
was  analyzed  in  [2],  where  it  was  shown  that  no  QR  algorithm  could  exist.  A  simpler  proof 
arises  from  noting  that  two  steps  of  LLT  is  equivalent  to  one  step  of  QR  in  the  positive 
definite  case,  and  so  the  question  is  whether  the  sparsity  pattern  of  To  =  LLT  is  the  same 
as  that  of  Xi  =  LT L\  this  is  easily  seen  to  include  only  tridiagonal  To. 
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